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I. INTRODUCTION 



Effective field theories have proved to be a powerful tool in a variety of applications, 
ranging over areas as diverse as nuclear physics 1], lattice QCD heavy flavor physics 
0|, and Bose-Einstein condensation 0]. One of the first applications of the technique 
was to the bound state problem in Quantum Electrodynamics (QED). As the approach 
incorporates QED effects as perturbations to a nonrelativistic Schrodinger problem, it is 
known as nonrelativistic QED (NRQED). 

Effective field theories in general provide a way of treating physics that involves multiple 
scales. Specifically in the case of atomic physics, NRQED deals with the fact that three 
scales-the rest mass of the electron, m; the average three-momentum of an electron, ma; and 
the electron binding energy, ma 2 -all_play significant roles in radiative and recoil corrections. 
While the Bethe-Salpeter equation [Q] or three-dimensional variants of it 

Baa 

provide a 

consistent framework to carry out calculations of atomic properties, the implementation is 
sufficiently complicated so that while first-order calculations were carried out in the early 
1950's for a number of systems, the next order calculations were not completely evaluated 
for over 40 years. 

This situation has radically changed since the introduction of effective field theory tech- 
niques for QED bound state calculations. A set of calculations using various implementa- 
tions of NRQED that complete the evaluation of order a 4 Rydbergs energy shifts have been 
carried out over the last few years: we note progress in helium [lQ], positronium [ijj, and 
muonium Q|. In addition to this work on energy levels, the order a 2 corrections to the 
decay rate of orthopositronium and parapositronium [3| have also been treated using 
NRQED. 

In this paper we will be concerned with the latter decay rate, r p „p s . The dominant decay 
mode of p-Ps is into two photons, and the associated decay rate was found by Wheeler 15j 
and Pirenne to be 



r^ Ps (2 7 ) = ]pY°^ = 27rcRooC * 3 = 8032.5028(l)jzs _1 . (1) 

The numerical value is determined using the 1998 CODATA adjustment of constants jl7|. 
with the 12 ppb uncertainty dominated by the uncertainty in the fine structure constant. 
In the following we will refer to this lowest order result simply as IV 
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The lowest order rate differs by 0.52 percent from the most accurate measurement 
which determines 



118], 



r;? Ps = 7990.9(1.7)^s- 1 . (2) 

The bulk of the difference is accounted for by the one-loop corrections to the decay rate, 
calculated by Harris and Brown which change the theoretical prediction by 0.59 percent 
to 

2 

r(lloop) = r {l + - 5)} = 7985.249 / us~ 1 . (3) 

The residual 0.07 percent discrepancy corresponds to a three standard deviation difference 
between experiment and theory at the one loop order, and is of a size compatible with order 
a 2 corrections. 

Corrections of this order arise both from two-loop corrections-the subject of this paper- 
and from the four photon decay of p-Ps. The experiment [3] measures the total decay rate 
of p-Ps, which in QED includes decays into any even number of photons. We note that exotic 
interactions could allow decay into an odd number of photons, and experiments looking for 
this kind of decay have put limits on the branching ratio arising from such interactions, 
specifically 2.8 x 10" 6 for p-Ps 37 Q and 2.7 x 10" 7 for p-Ps -> 57 Q The only effect 
that contributes at a non- negligible level is four-photon decay, which was first calculated in 
Ref. (2^. The highest accuracy determination of the rate, along with a calculation of first- 
order radiative corrections to it, is given in Ref. j^], where references to other calculations 
can be found. The result is 



r p _ Ps (4 7 ) = 0.274290(8)(^) 2 r {l - 14.5(6) Q}. (4) 



While this effect is well under the present level of experimental precision, we will include it 
in our final prediction. 

The remaining contribution from two and higher loops can be parameterized via 

r p _ Ps (2 + loop) = r f-2« 2 ln« + B 2 J-) 2 - ^ln 2 « + C— In a + d(-) 3 }. (5) 



The Q(a 2 ) logarithmic term was calculated in Ref. 24], and the leading 0(a 3 ) logarithm 



in Ref. |25J]. The coefficient of the subleading 0(a 3 ) logarithmic term 

C = — + 10 ln2 + — = 7.9189 (6) 
2 90 K ' 
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is also now well-established, with three different groups 2y, |27|, |28[ in agreement. The terms 
of order ct 3 ro are well below the experimental accuracy, but will be included in our final 
tally. The leading logarithmic term removes almost all the residual 0.07 percent discrepancy, 
so as long as the constant B 2l is not too large, theory and experiment are in agreement. The 
experimental error corresponds to a value of 39 for B 2l . Nevertheless, a direct calculation of 
this constant is desirable, both because it is possible that the constant is large, as is frequently 
the case in QED bound state calculations, and for comparison with future experiments of 
higher precision. As mentioned above, this calculation has been recently carried out using 
effective field theory techniques by Czarnecki, Melnikov, and Yelkhovsky jlJ ]. referred to 
here as CMY. After correction of one part of their calculation 2^1 , their result is 

5 27 = 5.1(3). (7) 

This gives for the total theory, including the four-photon decay, the prediction 

rf_ PS (CMY) = 7989.616(13)/is~ 1 , (8) 

which is in good agreement with experiment. The principal result of this paper is the 
confirmation of the previous calculation, but with higher precision, 

B 2l = 5.1243(33), (9) 

which leads to our main result, 

r^ Ps = 7989.6178(2)/is- 1 . (10) 

The plan of this paper is the following. In section II we explain our implementation of 
NRQED, which differs significantly from that of CMY. In section III we carry out a QED 
calculation of the scattering amplitude at threshold to two loops. In the concluding section 
we use this amplitude in a bound state calculation to determine B 2l , and we discuss the 
related o-Ps decay rate calculation. 



II. NONRELATIVISTIC QUANTUM ELECTRODYNAMICS 



The basic idea of NRQED is to take advantage of the fact that many of the complexities 
of QED are associated with the scale of the electron Compton wavelength, which in an 



atomic bound state is effectively a point interaction, so that in a bound state calculation 
one can account for most QED effects by introducing a delta function potential along with 
the usual relativistic perturbations. However, ultraviolet divergences are present when these 
interactions are treated in higher order, and the regularization of these divergences can be 
done in different ways. In addition, when carrying out scattering calculations, the infinite 
range of the Coulomb interaction leads to infrared divergences that can also be regularized 
in different ways. For these reasons, the details of an NRQED calculation can be quite 
different when done by different groups. We consider this to be an advantage, as agreement 
between different methods, such as will be found here, lends support to the reliability of 
these complex calculations. 

CMY handled both ultraviolet and infrared divergences with dimensional regularization. 
In the present work, we only use dimensional regularization to handle ultraviolet diver- 
gences in the QED scattering calculation of section III. We regulate infrared divergences by 
introducing a photon mass A = mX, and we cut off ultraviolet divergences in the NRQED 
calculation with a maximum loop three-momentum A. 

A detailed description of NRQED applied to the decay rate of orthopositronium can 
be found in Ref. [13]. A great deal of the analysis given there applies equally well to 
parapositronium, but for completeness we give here a complete, but somewhat abbreviated, 
discussion, referring the reader interested in more details to Ref. Q|. 

Because we regulate infrared infinities with a photon mass, it is necessary to set up 
a consistent set of interactions that incorporate this effect. These interactions are most 
conveniently expressed in the center-of-mass frame in momentum space, with an incoming 
electron momentum k and an outgoing electron momentum I. It is also useful to define the 
frequently occurring denominators 

D x (k)=k 2 + \ 2 (11) 

and 

D k = k 2 + 1 2 } (12) 

where 7 = ma/2. The NRQED interactions are depicted in Fig. 1. 

At lowest order, the NRQED Hamiltonian consists of the usual nonrelativistic kinetic 
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FIG. 1: The NRQED instantaneous potentials: (a) Coulomb, (b) relativistic mass increase, (c) 
Breit-Fermi, (d) four-fermion contact, (e) four-fermion derivative. 

energy together with a modified Coulomb potential, 

-> -» 47T0 1 

V c (k,l) = -—^ T -. (13) 
D x [k -I) 

The usual relativistic perturbations responsible for the fine structure of hydrogen are present, 
including the relativistic mass increase (RMI), 

V RM i(k,l) = -(2n) 3 5(k-l)— (14) 

and a term, which includes both relativistic corrections to the Coulomb potential and the ex- 
change of a transverse photon, denoted V B f- The full form of this "Breit-Fermi" interaction 
when a photon mass is present is 



-» -» Alio: 
V BF (k,l) = 



kxl\ 2 + \\ 2 \k + l\ 2 (k-l)x S--(k-l)x 



D\{k-l) D x {k-l) 



3.(kxl)-(S- + S+) 1 \k-l 



|2 



(15) 



2 D\{k — I ) ^D x (k-l) 
It includes the Darwin term and the spin-orbit interaction, though we note the latter inter- 
action does not contribute for S-states. 
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The next interaction we need for NRQED is the four-point interaction describing the 
decay process, 

V 4 (M) = V 4 °{l + - ei + (-) 2 e 2 } (16) 

where 

rrr 

The constants e\ and e 2 renormalize the interaction and will be determined from a matching 
calculation below. This interaction corresponds to a delta function in coordinate space and 
gives rise to an energy in which V4 is multiplied by the square of the wave function at the 
origin, m 3 a 3 /8n. The leading term is chosen so that the decay rate, given by T = —2Im(E), 
reproduces the known lowest order result. 

An important difference of the p-Ps calculation from the o-Ps calculation is the fact that 
the latter has a real contribution to V4 arising from one-photon annihilation, a channel not 
available in the present case. The last interaction needed for the calculation accounts for 
the fact that the annihilation is not exactly pointlike, which leads to a derivative term, 

v 4der (kJ) = ^t(P + P). (18) 



WM) = -^ 4 ° — — (19) 



We note the relationship 

2 P + P 

3 m 1 

happens to be identical to the behavior of the interaction for o-Ps arising from one-photon 
exchange, which contributes to the hfs of positronium; this fact will be used below. We now 
describe how these interactions are used in a bound state NRQED calculation. 



A. Bound state calculation 

We now apply standard Rayleigh-Schrodinger perturbation theory to second order. (The 
corresponding diagrams are shown in Fig. 2.) We are interested only in imaginary contri- 
butions to the energy, so one factor of V 4 or V ider must be present. The expression for first 
order perturbation theory in momentum space is 



where the wave function is 



/7 3 \V2 8tt7 






(d) 



FIG. 2: First and second order bound state perturbation theory contributions to the energy shift. 
They are the first order (a) contact and (b) derivative contributions, and the second order (c) 
relativistic mass increase and (d) Breit-Fermi contributions. Contributions (c) and (d) are shown 
separated into their zero-, one-, and many-potential parts. 

This gives the contributions 
from Fig. 2a and 



4L = — - (23) 



l(1) _ 2tAa 6 _ 1 ._ 7 
3tt C" 

from Fig. 2b, where A is the maximum value allowed for the magnitude of the three- 
momentum \p 2 \ or \f>i\. 

In second order we need to evaluate expressions of the form 



^ <0|V,|nxn|V,|0> 

Hlft — Uln 



which can be written in terms of a reduced Coulomb Green's function, 



The reduced Green's function in momentum space conveniently breaks into three parts, in 
which the electron either propagates freely, interacts once with a Coulomb potential, or 
interacts more than once. It is given by 



G(k,l 



m 



(2vr)V(fc-Z) — - 



m Ana m 



Dk D k \k-l\^D l 



-R(k,l 



(26) 



where 



R(k,l) 



327rm 7 3 r 5 4^ 2 4 7 



DID? 



with 



.4 



1 , a 2A ~ l 
+ - In A + - 

A 2 y/4A=T 



DM 



tan _1 \/4A - 1 



(27) 



(28) 



4 7 2 |A;-/| 2 

When the interaction is a relativistic mass increase (see Fig. 2c), second order perturba- 
tion theory gives 



p (2) , p (2) 



Aa t 



.ma' 



ln(- 



ma 

', 

32 



4tt 8 ' ^27- 

and when it is a Breit-Fermi interaction (see Fig. 2d), the contribution is 

. 7 ma 7 
I — 1. . 

'27' 



(29) 



E 



( 2 ) + Fj 
V 4 ,Vbf + ^ 



(2) 

V bf ,Va 



Aa 6 3ma 7 , / A 

-i i In ■ 

4tt 8 V 
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(30) 



where contributions of order ma 8 and higher are dropped. A table of integrals useful in 
arriving at these results can be found in Ref. 

The final result for the NRQED bound state calculation is then 



NRQED 



r a 
{l + - 



ei 



2A 
3m 



+ U h + ^ + 2- 2 ln(-)]}F , 



23tt 



A 



2 7 ^ 



(31) 



which has been converted to a decay rate in the usual manner. 



B. NRQED scattering calculation 

The next step in the analysis is to consider the imaginary part of the amplitude in NRQED 
for Bhabha scattering at threshold. In lowest order this is simply V®, with no contribution 
from V&der because of the presence of the momentum factors. Including the renormalization 
terms, we denote the overall effect of V4 as 

M i ° ) RQE o = V°{l + ^e 1 + Q 2 e 2 }. (32) 
9 



M}}> R qed(Vbf) = -V?{^ - ^) , (35) 



(a) (b) (c) (d) 



FIG. 3: The one-loop NRQED e e + — > e e + scattering graphs. 
In first order, a simple calculation gives the contributions 

M { A Q ed(V c ) = ^V 4 \^){i + % 1 }, (33) 

M« Q£D (W = ^4°(^-^), (34) 
-V°(— - — 
and 

M« Q ^(^ er ) = ^4°(-^ + ^p), (36) 
shown in Fig. 3 a-d, which sum to 

M^^ + ^-f^). (37) 

In the next section we will require this amplitude equal the one-loop QED amplitude to 
determine e\. To determine €2 a two-loop calculation is required. The calculation is almost 
identical to the scattering calculation carried out for hyperfine splitting in section 3 of Ref. 
[l3l |. where twelve diagrams (labeled a-1) were analyzed. The last two of them, k and I, 
involved a real contribution to V4 not present for p-Ps, so they are dropped. The other ten 
(see Fig. 4) are all proportional to the lowest order amplitude and can be taken over directly 
to p-Ps by simply using the lowest order amplitude for p-Ps instead of hfs. The only other 
change needed is the modification of the three diagrams which involve the Vbf interaction, 
which is different for the two states. However, the difference is simply 

87ra \k — 1 1 2 

SV BF = o - - , 38 

3m 2 D x (k-l) 

which leads to the modifications 

V=©V(^"^). (39) 

10 




FIG. 4: The two-loop NRQED e e + — > e e + scattering graphs. 

and 

(«) 

Finally, graphs z and j are derivative terms, but as mentioned earlier, the ratio of the 
derivative term to the lowest order matrix element is the same for hfs as it is for p-Ps, 
S0 the res* frora Ret. Q caa be ta.ee ove r directl, The Beat result for tee two-loop 
contribution to threshold scattering is 

,A2) /a\2 Trn /7r 2 (21n2 + 1) Auk n 9l /An 47r 2 x . . 

M^ D =(-) ^-^ + 2^^) + —), (42) 

and the NRQED amplitude that is to be compared to the QED scattering calculation is 

M NRQED = MtfRQED + M NRQED + M NRQED- ( 4 3) 
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III. QED SCATTERING CALCULATION 



We are interested in the imaginary part of the scattering amplitude for an electron and 
positron to annihilate into two photons, and the two photons then to create an electron 
and positron. While this amplitude at threshold is simply a constant, in order to determine 
the derivative term in NRQED we also need the behavior slightly above threshold. If we 
assign the momentum k to the incoming electron and I to the outgoing electron, with the 
positrons having opposite momentum, a straightforward calculation leads to the above- 
threshold matrix element 

*&(*.0 = v?(i-ft£). («) 

This requires that NRQED have the forms for V4 and V^er used in the previous section. At 
threshold we have simply 

<L = V 4 °. (45) 

We will follow the convention of listing results for a given n-loop contribution to the scat- 
tering amplitude 

<L = ©>4°' (n) (46) 
by giving / instead of the full amplitude. So for the zero-loop contribution we have 

J(°) = 1. (47) 

Continuing to the one-loop calculation, there are three diagrams that must be considered, 
which we call the self-energy (SE), vertex (V), and ladder (L). They are shown in Fig. 5. 
The self-energy and vertex results are 

/ 5E = 2mA + 41n2 + l + 0(A), (48) 

and 

l v = -41nA + — -41n2-4 + 0(A). (49) 

The infrared logarithms present in I$e and Iy are artifacts of the renormalization process. 
The ladder diagram has an infrared singularity associated with binding that we denote as 
Ib, defined by 

/ B = ^ + lnA-l-^ + 0(A 2 ). (50) 
A o 
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(SE) 



(V) 



(L) 



FIG. 5: One-loop QED graphs contributing to the parapositronium decay rate. They are (SE) 
the self-energy graph, (V) the vertex graph, and (L) the ladder graph. 

The ladder diagram turns out to be proportional to Ib, which will also prove a useful factor 
in the two-loop calculation. We find 



I L = 2I B + 0{\). 

A simple summation of the one-loop contributions gives the overall result 

/ (1) = f +A(X), 



with 



7T 



A(X) = — -5 + EX. 



(51) 



(52) 



(53) 



Here EX represents the uncalculated order A contribution. The value of E does not affect 
the final result. 

The two-loop calculation involves eight sets of graphs. In this part of the calculation 
ultraviolet divergences are regulated by working in n = 4 — 2e dimensions. When doing so a 
common factor of [(47r / u 2 /m 2 )e~ 7E ] 2e is always present (where /i is a mass scale introduced 
in the process of dimensional regularization, and 7# is the Euler gamma constant), but we 
leave it unwritten. It can be taken to unity after renormalization. 

We begin the discussion with two-loop corrections to a photon vertex, shown in Fig. 6a. 
In diagrams a2 and a5 the self-energy subdiagram is understood to be accompanied by a 
self-mass counterterm. This convention differs from CMY, but we have been able to show 
agreement with their results if the counterterm is included. The unrenormalized results are 
given in Table I. The total class a unrenormalized result is 
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FIG. 6: Contribution to the p-Ps decay rate from groups a (the two-loop vertex corrections), b 
(the two-loop self-energy corrections), and c (products of one-loop vertex and self-energy parts). 

I' a = ± + -J^ - In 2 - 1) + 2.8841255(52). (54) 
This is renormalized according to 

I a = I' a -i 1 I' v -A{i 2 -i\)I^\ (55) 

where 

1 2 t 2 

/(, = - + - 4 In 2) + e(^C(3) + y + 4 In 2 2 - 8 In 2) + 0(e 2 ) (56) 
is the unrenormalized one-loop vertex correction (see Fig. 5); L\ = (a/7r)£i, where 

h = 1 + (in A + l) + e(- In 2 A + ^ + 2) + 0(e 2 ), (57) 

is the one-loop vertex renormalization constant; and L 2 = (a/ir) 2 ^, where 

- 1 1/lnA 13x /l lSl , , 3 W . 2l 157tt 2 49 x . , „. 

= 32? + + S) + (j A + A + 2 C(3 > " * ln 2 + "ToF - l2i) + < 58 > 

is the two-loop renormalization constant (with the vacuum polarization effect excluded). 

Because we treat vacuum polarization separately, it is important to again emphasize that 

the renormalization constant above also does not include the effect of vacuum polarization. 

Our result for the two-loop renormalized class a contribution is 

I a = -2 In 2 A - lnAJy(A) - 1.966447(6), (59) 
14 



where Iy(A) is the one- loop vertex factor of Eq. 

Turning to the two-loop self-energy diagrams, shown in Fig. 6b, we again present results 
with self-mass counterterms understood. The unrenormalized results are given in Table II 
and are in agreement, after accounting for the different convention, with CMY. The total 
class b unrenormalized result is 

K = + 7(^2 - A) _ 1.6727038(15). (60) 

The corresponding renormalized expression is 

h = I' h ~ ill'sE + 2(4 - ti) 

= In 2 A - In Mse(X) - 2.0277430(15), (61) 



I'se = ~y + (4 In 2 - l) + - 4 In 2 2 + 6 In 2 - 2) + 0(e 2 ), (62) 



where I' SE is the unrenormalized one-loop self-energy correction 

r 2 

2e ' V ' "V 24 

and Ise{^) is the renormalized one-loop self-energy factor of Eq. 

The diagrams with two separate one-loop corrections shown in Fig. 6c are simple to 
evaluate and can be done entirely analytically. The breakdown of the calculation is given 
in Table III and agrees exactly with CMY. Interestingly, the total class c contribution is 
unchanged by renormalization. It is 

I c = — + ivr 2 In 2 - — + ln 2 2 - In 2 + -. (63) 
c 128 4 8 4 K J 

The least complicated contributions are those, shown in Fig. 7d, which have no ultraviolet 
or infrared divergences. The results, given in Table IV, can be directly compared with the 
diagrams called D\± and D% by CMY and are in agreement. For the class d total we find 

I d = -0.87783(7). (64) 



The vacuum polarization diagrams of Fig. 7e were evaluated in Refs. |30l |3jj , with result 

I e = 0.4473430(6). (65) 

As with o-Ps, by far the most difficult diagrams to treat were those of class f, with at 
least one binding photon exchanged. As before, the most infrared singular contributions are 
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FIG. 7: Contribution to the p-Ps decay rate from groups d (infrared and ultraviolet finite contri- 
butions), e (vacuum polarization corrected one- loop graphs), and f (ladder graphs). 

the ladder and crossed- ladder fl and f2. They were treated just as in the o-Ps calculation 
by subtracting off the most infrared singular part of each (called flaO and f2a0), which were 
evaluated as a unit. We found that 

Ifi = Ifiao + Ib{X)Ils(X) - 3tt 2 In A - 55.2277(15), (66) 
I}2 = / /2 ao + 7r 2 hiA + 30.8346(7), (67) 

where (see Eq. (jUJ)) Ils = II — 21 b = 0(A). The leading binding singularity is contained 
in 

I/iao + //2ao = 2/ a0 = + 2(l B - + In 2 A - 21nA + 2A , (68) 

where 

A = -3.2036(13) (69) 

as given in Eqs. (3.37) and (3.86) of Diagrams f3 and f4 were evaluated using Yennie 
gauge for the vertex or self-energy photon and Feynman gauge for the ladder photon. This 
procedure is consistent since the sum of the two is independent of the gauge of the vertex or 
self-energy photon. An all-Feynman-gauge evaluation of f3 and f4 proved to be numerically 
difficult. The class f results are given in Table V. The total class f contribution is 

I f = 27r ^ n2 + ^A(X) + In 2 A + (iy(A) + Ise(X) - 2vr 2 ) In A - 35.0357(32). (70) 
A A 
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Class g is the set of energy shift contributions having one-loop corrections (SE, V, or L) 
both before and after the annihilation into two photons. This contribution is easily evaluated 
as the square of half the one-loop correction: 



1k» 

2 



2 7T 2 I ,,,, /7T 2 5\2 



Finally, class h represents the effect of light-by-light scattering on the annihilation pho- 
tons. This contribution was evaluated in Ref. [29J to be 

I h = 1.29392(4). (72) 

The total two-loop contribution is then 

M 7T 2 (21n2 + 1) 2tt 9l a 

1 = x2 + -r^(A) - 2tt 2 In A + B 2 , (73) 

A A 

with 

B 2 = -35.2881(33). (74) 

The total QED scattering amplitude at threshold through two-loop order is the sum of the 
Mq ED for i equal to 0, 1, and 2. 

Most of the two-loop results can be cross-checked against corresponding results of CMY 
as given in Ref. Q] . This comparison is detailed in Table VII. The a graphs can be compared 
directly to the corresponding CMY graphs except for a2 and a5, whose CMY counterparts 
Du and D 12 require the inclusion of the one- loop self- mass 5m^ as given in CMY Eq. (48). 
(For ease of notation we take m — > 1 in the CMY expressions for 8m.) The CMY coun- 
terparts of the b graphs require one- and two-loop self-mass insertions. It is necessary to 
isolate the "crossed rainbow" and "double rainbow" contributions to 5rrS 2 \ which are 

SmPn = 4-!:-12C(3) + 8vr 2 ln2-^ + i (75) 



as given in Ref. 
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*»E. = £ + 5 + £ + X- (76) 

f. 3^| but adapted to the conventions of CMY. The c graphs were done 
analytically by both groups, and agree exactly once appropriate self-mass insertions are 
made. The d graphs are directly comparable. The e (vacuum polarization) graphs were 
taken from Ref. 3l( by both groups. The f graphs, which involve the binding interaction, 
were the most difficult to compare. In fact, the double ladder and crossed ladder graphs 

17 



fl and f2 could not be checked against their CMY counterparts D4 and D\ because the 
methods for regulating the binding singularities were completely different. We evaluated 
graphs f3 and f4 with Yennie gauge vertex and self-energy photons, but the sum is directly 
comparable with D 2 + D 5 of CMY after the appropriate self- mass correction is made to their 
result. The comparison for f5 and f6 is a bit more delicate. For the sum f5+f6 we found 

I fs + J /6 = (/y + I se )Ib + 1.29396(4), (77) 

where ly and Ise are the renormalized vertex and self-energy corrections. Now CMY cal- 
culate unrenormalized graphs, so we revert to the unrenormalized contribution as well: 
ly + Ise I'v + I'se-i where in CMY notation 

I' v + I' SE = 4(S 2 + S 3 + Sm^B 1 ). (78) 

Now 1b is half the contribution of the one-loop ladder graph, which in CMY terms is 
Ib = (4Si)/2 = 2 Si. Note that we must regulate the infrared divergence via dimensional 
regulation for purposes of comparison, since that is the method of CMY. Now the sum 
J/5 + J/6, with the substitutions for ly + Ise and Ib discussed above, is given by 

J /5 + //6 -L + - l) + 5.84851(4), (79) 

while the corresponding CMY contribution is 

D 6 + D 7 + 5m {1) C 7 = + - - l) + 5.850(12). (80) 



IV. RESULTS AND CONCLUSION 



The NRQED couplings e± and e 2 are determined by requiring agreement between the 
imaginary parts of the scattering amplitude as calculated by NRQED and QED. The 
NRQED version of this amplitude, given in Eq. PHj) . is 



M NRQED — ^4°{ 

+ (-) 2 



a 

1 + -e x 

7T 



«\2 



+ (") 

7r 2 (21n2 + 1 



e 2 



+ 



a 



71 



2tt 



A v 



(l 



a 



47rA 



7T 

A 



-d 



2A 7ttA 
3^ + ~12 



^ + 27r 2 ln(^) + — }. 



A 2 



3A 



2A' 



The corresponding QED version is 
Mqed = V°{l + 
+ 



cv 



7T 



2tt 

T 



+ A(X) 



ax2r7r 2 (21n2 + 1) 

IT' - 



+ ^A(A)-27r 2 lnA + J B 2 ]}. 



A 2 



A 



(81) 



(82) 
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Matching at the one-loop (order a) level yields 

2A 7ttA 

To order aT , this value of e± used in Eq. (|3"Tj) cancels the ultraviolet divergent term and 
reproduces the Harris-Brown result ^| after taking the limit A —>■ 0. Matching at the 
two- loop (order a 2 ) level gives 

e2 = - 2 * 2 "<^)-Y + ^ (84) 

Care is required with the photon mass A in this evaluation. While it is ultimately taken to 
vanish, we keep it in e\ because that factor enters as a factor of 1/A terms. However, while 
e 2 also has terms proportional to the photon mass, they can be dropped because we do not 
need the next renormalization constant e%, which enters in higher order. We note that the 
uncalculated 0(A) part of A(X) cancels from the calculation. 

With ei and e 2 now determined from the matching calculation, we return to the bound 
state NRQED calculation and use these values in Eq. (|31j) . Taking the limit A — > then 
gives our final result for the two photon decay rate through two-loop order: 

2 

Tnrqed = {l + - 5) + (^) 2 [-2vr 2 ln« + B 2l ]}v , (85) 

where 

65n 2 

# 27 = # 2 + __ + 2tt 2 In 2 = 5. 1243(33). (86) 

As noted in the introduction, this is consistent with the CMY result. Inclusion of all known 
contributions leads to 

r p _ Ps = \l + A- -2a 2 lna + B(-) In 2 a + C—lna + D{-) \T , 87 

L 7T V 7T / 27T 7T V 7T /J 

where A = tt 2 /A — 5, 

B = B 2l + B^ = 5.3986(33), (88) 
C = 7.9189, and D is uncalculated. Our final numerical result is 

r p _ Ps = 7989.6178(2) / us~ 1 . (89) 

The uncalculated D term makes a contribution of O.OOOlOZ^/xs -1 . Numerical contributions 
of the various terms in T p _p s are detailed in Table VIII. While the experimental precision is 
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relatively high for a decay measurement, it is far too low to be sensitive to the precise value 
of B, only disfavoring a very large value on the order of 50 or more. The very short lifetime 
of the state makes it unlikely that significant improvements in the experimental precision 
can be achieved in the near future. 

We consider one of the most important aspects of the present calculation to be its impact 
on the decay rate of o-Ps. We have applied exactly the same methods to both decay rates. In 
the o-Ps case, while no independent check of our calculation was then, or is now, available, we 
did use the same implementation of NRQED to calculate a contribution to the ground state 
hyperfine splitting of positronium that had been determined using Coulomb gauge Bethe- 

1 1 n 

Salpeter methods |33J] and additionally an independent NRQED approach [34], finding good 
agreement. This indirect confirmation of the validity of our methods has now been further 
buttressed by the fact that we have agreement for the decay rate of p-Ps with CMY, who, as 
we have emphasized above, use an entirely different implementation of NRQED. Therefore, 
the confrontation of QED with experiment in this system, summarized as 

r* Ps [13] = 7.039979(ll) / us~ 1 

r^ p Ps [35] = 7.0514(14) / u^ 1 

T e % B [36] = 7.0482(16)/is~ 1 

r^ p Ps [37] = 7.0398(29)/is" 1 (90) 

remains in an unresolved state, and strongly indicates the need for further experimental 
work. 
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TABLE I: Unrenormalized contributions to the parapositronium decay rate from class a. Class a 
consists of two- loop vertex corrections. 



diagram 


1 




1 

e 




1 


al 


1 

8 


7T 2 

16 


- In 2 + 


3 
16 


-0.3237469(14) 


a2 


_ 1 

8 


16 


+ ln2 - 


. J_ 
16 


-0.3245758(11) 


a3 


1 

8 


7T 2 

16 


-In 2 + 


5 

16 


1.4769869(20) 


a4 







l 

2 




3.492237(3) 


a5 


1 

8 


7T 2 

16 


+ ln2 - 


3 
16 


-1.330000(3) 


a6 


1 

8 


7T 2 

16 


-In 2 + 


3 
16 


-0.1067757(12) 


total 


1 

8 


7T 2 

16 


- In 2 - 


3 
16 


2.8841255(52) 



TABLE II: Unrenormalized contributions to the parapositronium decay rate from class b. Class b 
consists of two-loop self-energy corrections. 



diagram 


l 

-? 


1 

e 


i 


bl 

b2 


l 

8 

1 

16 


21n2-£ 
-ln2 + £ 


-1.4022755(10) 
-0.2704283(10) 


total 


1 

16 


ln2 -i 


-1.6727038(15) 



TABLE III: Unrenormalized contributions to the parapositronium decay rate from class c. Class 
c consists of products of two one-loop vertex or self-energy corrections. 



diagram £ ^ ^ \ it 4 C(3)vr 2 ln2 it 2 In 2 2 In 2 1 



cl 


l l 

4 16 
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7 
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1 

24 
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-1 
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1 1 

8 16 


-1 





1 

128 


7 

8 


1 
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1 

16 


2 


-2 





c3 


1 o 


-2 


1 

2 











7 

48 


6 


-5 


5 
4 



total OOOOj^gO \ -\ 1 -1 \ 
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TABLE IV: Contributions to the parapositronium decay rate from class d. Class d consists of 
ultraviolet and infrared finite two-loop corrections. 



diagram 


1 


dl 


-0.80394(3) 


d2 


-0.07389(6) 


total 


-0.87783(7) 



TABLE V: Contributions to the parapositronium decay rate from class f. Class f contains the 
ladder graphs. 



diagram 


Ib(X) 


vr 2 In A 


1 


fl-flaO 


IlsW 


-3 


-55.2277(15) 


f2-f2a0 





1 


30.8346(7) 


f3 








-10.05855(11) 


f4 








3.9965673(13) 


f5 


MA) 





1.94968(4) 


f6 


Ise(X) 





-0.655722(5) 


total 


Ils(X) + I V (\) + Ise(X) 


-2 


-29.1611(17) 



TABLE VI: Renormalized two-loop QED contributions to the parapositronium decay rate by class. 



class 


7T 

A 


In 2 A 


In A 
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a 





-2 


-MA) 


-1.966447(6) 


b 





1 


-Ise(X) 


-2.0277430(15) 


c 











1.274886 
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-0.87783(7) 


e 











0.4473430(6) 


f 2 In 2 


A(X) 


1 


Iv(X)+Ise(X)-27t 2 


-35.0357(32) 


g 1 


A(X) 








1.603514 


h 











1.29392(4) 


total 2 In 2 + 1 2A(\) 





-27T 2 


-35.2881(33) 
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TABLE VII: Cross-check of results between this work (AMFS) and CMY. The 1/e 2 and 1/e terms 
agree analytically in all cases, only the finite parts are displayed. The comparison for f5+f6 is 
subtle and is discussed in the text. 



AMFS 


CMY 


AMFS 


CMY 


i -\ ■ f i" 

Difference 


Graph 


Graph 


Result 


Result 


AMFS-CMY 


al 


#8 


-0.3237469(14) 


-0.324(1) 


0.000(1) 


a2 


Dn + (W^Cn 


-0.3245758(11) 


-0.33(3) 


0.01(3) 


a3 


A) 


1.4769869(20) 


1.475(50) 


0.00(5) 


a4 


£>io 


3.492237(3) 


3.488(2) 


0.004(2) 


a5 


Z?i2 + Sm^C 12 


-1.330000(3) 


-1.27(12) 


-0.06(12) 


a6 


D13 


-0.1067757(12) 


-0.12(1) 


0.01(1) 


bl 


Die + ^ci^i 


-1.4022755(10) 


-1.403(4) 


0.001(4) 


b2 


L>i7 + Jm^^Cir + 5m K D ' R B 1 


-0.2704283(10) 


-0.271(1) 


0.001(1) 


cl 
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analytic 
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analytic 


analytic 





c3 


D 15 + 5m^C 15 + (5m^) 2 B 2 


analytic 
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dl 


D u 


-0.80394(3) 


-0.804(2) 


0.000(2) 


d2 


D 3 


-0.07389(6) 


-0.074(2) 


0.000(2) 


f3+f4 


D 2 + D 5 + 5mV>C 5 


-6.06198(11) 


-6.03(20) 


-0.03(20) 


f5+f6 


D 6 + D 7 + 5m^C 7 


5.84851(4) 


5.850(12) 


-0.001(12) 
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TABLE VIII: Numerical values of contributions to the parapositronium decay rate. 

term contribution (in /is -1 ) 

1 8032.5028(1) 

-47.2534 
-2a 2 In a 4.2092 
5(f) 2 0.2340(1) 
-%^ln 2 a -0.0361 

Z7T 

C^lna -0.0387 

7T 

L>(f ) 3 0.00010 D 

total 7989.6178(2) + 0.00010D 
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